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Abstract 

We present a new first-principles formalism for calculating forces for optically excited electronic states 
using the interacting Green's function approach with the GW-Bethe Salpeter Equation method. This 
advance allows for efficient computation of gradients of the excited-state Born-Oppenheimer energy, al- 
lowing for the study of relaxation, molecular dynamics, and photoluminescence of excited states. The 
approach is tested on photoexcited carbon dioxide and ammonia molecules, and the calculations accu- 
rately describe the excitation energies and photoinduced structural deformations. 

Calculation of optically excited electronic states and spectra has recently become possible through use of 
the interacting two-particle Green's function within the first-principles GW-Bethe Salpeter Equation (GW- 
BSE) formalism [1, 2], making possible the study of photoinduced structural change and photoluminescence. 
However, within this methodology, one calculates optical properties at fixed ionic positions: barring inefficient 
finite-difference schemes or shrewd guesses, the direction in the multidimensional space of ionic configurations 
best optimizing the geometry is unknown. The problem is of practical importance since structural changes 
due to optical excitation are general phenomena which cause atomic rearrangement or dissociation and 
change the structure or symmetry of defects. The possibility of efficient ab initio calculation of excited-state 
forces opens new doors to reliable modeling and study of such changes. 

To this end, we present a new formalism for calculating excited-state forces within the GW-BSE approach. 
We develop the theory for the force calculations and the approximations requisite to render computations 
tractable and then carry out practical tests for two molecules. We also compare our results to those of 
constrained density functional theory (CDFT) [5]. Our development of force calculations parallels recent 
quantum chemistry advances (e.g. [3]) where analytical forces can be calculated for excited-state methods 
(CIS, CIS(D), EOM-CCSD). These methods provide an overall accuracy similar to the GW-BSE method [3, 
4]. However, the GW-BSE method scales as TV 4 (N is the number of atoms in the system) whereas these 
methods scale as iV 6 or worse. In addition, quantum chemistry methods are much more difficult to apply to 
the bulk properties of solids whereas the GW-BSE method scales equally well in the bulk limit. 

Within the ab initio GW-BSE approach, the ground-state electron density, total energy, forces, and 
single-particle states \i) d ^ t and eigenvalues s/ 1 are obtained using density functional theory (DFT). We then 
construct the RPA dielectric function e _1 and the screened interaction W = e~ 1 v c , where v c is the Coulomb 
interaction v c {f, f 1 ) — l/\f—r'\. Quasiparticle excitations are found by using the GW approximation to the 
self-energy, E ~ iGW [6]. The effective quasiparticle Hamiltonian H qp is 

H«p = T +V sc + (±-V xc ), (1) 

where T is the kinetic operator and V sc = Vi on + Vh + V xc is the sum of the ionic, Hartree, and DFT mean- 
field exchange-correlation potentials. We solve the Dyson equation H qp \i) = Ei\i) to obtain quasiparticle 
energies £i and eigenstates \i). 

Two-particle excited-state properties are obtained by solving the Bethe-Salpeter equation of the two- 
particle Green's function. We employ the "standard" positive-frequency version of the BSE eigenvalue 



equation and restrict to static screening [1, 2]. The BSE eigenvalue equation is 



y^, H cv S dv' A c'v' - ^SA S CV1 (2) 
c'v' 

with Hf v S ,tf v < = (e c - £v)S cc 'Svv' + K cv . dv , . (3) 

Here c labels unoccupied (conduction) states and v labels occupied (valence) states, Af v is the electron-hole 
amplitude for a quasihole in state v and a quasielectron in state c, S labels an excited state, fig is the 
excitation energy, and H BSE is the effective electron-hole Hamiltonian. The excited state energy is given by 
Es — Eo + where E is the ground state energy. K is the electron- hole interaction kernel with matrix 
elements 

K cv>c , v ,= y c (l)*t;(2)S(1234)c'(3)t; / (4)*d(1234), (4) 

where — The kernel S is 

5(1234) = -5(13)5(24)W(12) + 5(12)5(34)« c (13) . (5) 

The static screening approximation means that W = e~ 1 v c uses the static dielectric function e _1 (u; = 0) in 
Eq. (5). We note that the frequency dependent e _1 (w) is used in the GW quasiparticle calculations. We use 
the standard normalization J2 C v l^fJ 2 = 1- 

Our aim is to compute excited-state forces, i.e. derivatives of Es versus the 37V ionic coordinates R, 
denoted as OrEs- The derivatives have two parts 

OrE s = d R E + d R n s ■ 

DFT provides the ground-state derivatives 8rE [8]. The BSE provides expressions for fig, and we compute 
OrQs directly. Using Eq. (3) and the normalization condition for A^ v , we have 

d R il s = A s cv *A s c , v ,d R H*%, (6) 

cv,c'v' 

where 

d R Hc V s c ? v i = {dRe c — dRe v )6 cc '5 vv > + dRK CVtC i v > . (7) 

The derivative 8rH bse contains two types of terms, those involving &R£i and those involving d R K. Below, 
we adopt two physical approximations to render computations tractable. 

Since H qp \i) = Ei\i) is a standard eigenvalue equation, standard first order perturbation theory yields 



d R£l = (i\d R Hip\i) 








pn ee (j| W)} = j u\e K Api) if ; I • (8) 

Evaluating dRH qp = 8rV sc + 9^(S — V xc ) is burdensome as OrT, contains the derivatives OrG and Ore^ 1 , 
both prohibitive to calculate. We choose a physical approximation: as discussed above, V sc contains the 
dominant ionic, Hartree, and mean-field exchange-correlation potentials which bind solids and molecules. 
The weaker correction E — V xc consists largely of a constant "scissors-shift" with a weak dependence on R. 
Therefore, we approximate dR,H qp = OrVsc- To calculate 8rV sc , we employ density functional perturbation 
theory [7]: an auxiliary quadratic functional is minimized, and at its minimum OrVsc is easily found. We 
perform 37V minimizations for the 37V choices of R. 

Considering 8rK and examining Eqs. (4) and (5), there are two distinct types of derivatives: those 
containing 8r\i) and those containing 8rW. For the first set, we employ Eq. (8) and sum over intermediate 
states to reach convergence. For the second set, we now argue (and numerically verify below) that they are 
negligible. Specifically, in the GW-BSE formalism, W = e~ x v c where e = I — v c P and P is the polarizability. 
Within the RPA approximation that we employ, P = iGG is a function of G alone and thus so are e and 
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W, so the chain rule yields 8rW(12) — J jjj^ • 9flG(34)d(34). In deriving of the interaction kernel 5, we 
assume that 5W/8G = [2]. Judging from the success and accuracy of the BSE based on this assumption, 
we conclude that we may set 8rW = 0. Thus, we arrive at the following expression for OrK: 

3 

P R , K ■ i 4- P R ,* K (9) 

Taken together, Eqs. (6-9) provide us with an explicit expression for Or£Is and hence OrEs- 

We now present test applications to verify the accuracy of our approach for molecules for which high 
quality excited-state observations are available. We begin by considering the first singlet excited state of 
carbon monoxide (CO) which has one degree of freedom R. 

We carry out density functional calculations using the plane-wave pscudopotcntial method within the 
local density approximation (LDA) [8]. We use Kleinmann-Bylandcr pscudopotcntials with s and p projectors 
(p local) for C and O with cutoff radii of r c = 1.3 a.u. [9] and expand the electronic states to a plane wave 
cutoff of 70 Ry. Our periodic supercell is a 7 x 7 x 7 A 3 cube. We sample the Brillouin zone at k = 0. We 
perform GW calculations using the generalized plasmon-polc model [6]. We include 600 bands in the GW 
calculations to converge absolute quasiparticle and ionization energies. In the GW-BSE computations, we 
truncate the Coulomb interaction beyond 3.5 A to avoid spurious periodic image interactions [2]. With these 
parameters, all reported energies are converged to 0.05 eV. For comparison, we also calculate excited-state 
properties with the constrained-LDA (CLDA) method: we occupy the LUMO with an electron taken from 
the HOMO and iterate to self-consistency. 

We mention two technical issues: (i) off-diagonal matrix elements of H qp in the DFT basis are required: 
their neglect changes fig by ~ ±0.2eV, as noted previously [10]; (ii) in the DFT supercell calculations, the 
Coulomb interaction has infinite range: this creates an ambiguity in the vacuum level and shifts quasiparticle 
energies by a constant. By varying the supercell volume V, we find the shift proportional to 1/V and 
extrapolate to V — > oo. 

Fig. 1 displays the energy of the ground state and first excited state of CO as a function of the bond 
length. Table 1 lists calculated and experimental properties of the ground state (A 1 E + ) and excited state 
(A 1 !!): equilibrium bond lengths R e , harmonic vibration frequency ui e , ionization energy IP (the HOMO 
energy for the LDA), and the minimum-to- minimum transition energy T e . Results are also presented for 
representative quantum chemical methods for which forces have been developed. The ground-state LDA 
R e and u e are in good agreement with experimental results, and the quantum chemical results are slightly 
superior due to the better treatment of correlation. While the LDA HOMO energy lies far from the ionization 
energy, the GW results remove this error. The BSE results for the transition energy T e to the excited state 
agree with experiment with an error typical of the method [2] whereas the CLDA transition energy is off by 
1 eV. For R e and uj e both the CLDA and BSE perform equally well; while we do not have an explanation for 
the size of the error, the fact that for this particular excitation the CLDA produces an u e of similar quality 
to the BSE one is explained by the observation (see below) that both methods produce similar variation of 
force versus bond length. The BSE results for T e are of equal or better quality than the quantum chemical 
results whereas for R e and co e they are slightly worse. 

Turning to the forces, we wish to know how accurately we can compute OrEs, i.e. the slope of the curves 
in Fig. 1. Fig. 2 presents the calculated forces, as formulated above, along with "exact" forces obtained from 
two-point finite differences of the BSE energies. We also compute forces by assuming 8rK = in Eq. (7), a 
quasiparticle-only treatment and the closest analogue to the single-particle CLDA. Both these single-particle 
methods predict equal and opposite forces for the C and O atoms, signaling that dRH qp = dRV sc is a good 
approximation. Interestingly, both produce similar forces that depart from the "exact" forces by essentially 
a constant. Hence, both methods predict well the variation of the force versus R, an a posteriori verification 
of our physical intuition that changes in the mean-field potential V sc dominate. 

When we include contributions from 8rK , the forces improve markedly. Since we assume that 8rW = 0, 
the forces on C and O are no longer exactly equal and opposite. However, as Fig. 2 shows, the deviations are 
quite small on the relevant scale: as remarked above, the approximation 8rW = is excellent in practice. 
If we subtract the unphysical net force on the center of mass (i.e. averaging the C and O values in Fig. 2), 



3 



the calculated force become essentially "exact" . 

CO has a single degree of freedom allowing for careful study. However, calculation of forces is truly useful 
when there are many degrees of freedom and one does not know, a priori, which are the relevant ones for 
a given excited state. We now consider the first singlet excited state of ammonia, NH3, a molecule with 
six degrees of freedom. We employ identical methods as in the case of CO and provide key parameters: s 
projector for H (r c =0.8 a.u.); s and p projectors for N (p local, r c =1.0 a.u.); identical supercell, k-points, 
and Coulomb truncation radius; a 50 Ry cutoff for the wave functions; 600 bands in the GW portion; all 
energies converged to 0.05 eV. 

Table 2 lists properties of the LDA ground state. Starting with this ground state, for the BSE we 
consider the first singlet state, and for the CLDA we promote an electron from HOMO to LUMO. We 
compute excited state forces and perform relaxations until bond lengths are converged to 0.01 A, bond 
angles to 1°, and transition energies to 0.01 eV. Table 2 presents results for the relaxed excited state using 
CLDA, GW-BSE, and representative quantum chemical methods. The BSE excitation energy compares well 
with experimental and quantum chemical values. (The flattening of NH3 is along the famous "umbrella" 
mode.) 

Based on these two cases, we have verified the accuracy of the BSE method for calculation of excited state 
energies and geometries. Our approach provides excited-state forces which are, to an excellent approximation, 
the derivatives of the BSE energies. Intriguingly, for these two cases, CDFT yields inferior excitation energies 
but predicts geometries of comparable quality to the BSE ones. While encouraging, there are a number of 
serious problems with wider use of CDFT. Use of CDFT is straightforward when the excited state is composed 
mainly of a single configuration: in the cases above, the HOMO-LUMO combination has probability above 
90% for the excited state. However, for higher excited states or larger systems, we can have multiple 
configurations, something we can not know without solving the BSE. Therefore, we believe that CDFT 
may be a useful guide in certain circumstances, but results thus obtained must be carefully tested by more 
sophisticated methods. 

In brief, we present an ab initio formalism for calculating excited-state forces within the GW-BSE method 
as well as approximations allowing for computational tractability. We compute the photoexcited properties 
of molecules and verify the accuracy of (a) the GW-BSE formalism for describing the excited-state energies 
and structural relaxations, and (b) the forces as per our formalism. The calculations are as accurate as 
leading methods used in quantum chemistry (for which analytical force calculations are available) while 
scaling significantly better with system size and being easily applicable to the bulk. 
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Figure 1: ground-state LDA (circles) and A 1 !! first excited-state energies (triangles are GW-BSE, 

stars are CLDA) versus the bond length R for CO. The continuous curves are polynomial fits. 
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Figure 2: Absolute magnitudes of A 1 !! excited state forces for CO. Circles are "exact" GW-BSE forces. 
Crosses and pluses are GW-BSE forces on C and O. Triangles are GW-BSE forces with 8rK = 0. Stars are 
CLDA forces. The dashed curve is a guide for the "exact" forces. 
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LDA 

GW 
MP 2 [3] 
CCSD [3] 
Expt. [11] 



CLDA 
GW-BSE 
CIS [3] 
EOM-CCSD [3] 
Expt. [11] 



Ground state 

R e (A) u e (cm- 1 ) IP (eV) 

1.13 2050 9.1 

14.1 

1.133 2151 
1.124 2243 

1.128 2170 14.01 

Excited state (A 1 !!) 

R e (A) u e (cm- 1 ) T e (eV) 

1.21 1720 7.02 

1.26 1290 8.32 

1.21 1633 8.83 

1.22 1593 7.91 
1.24 1518 8.07 



Table 1: Ground state and excited state data for CO: equilibrium bond length R ei harmonic vibrational 
frequency u e , ionization potential IP, and — » A 1 II minimum-to- minimum transition energy T e . 



Ground state {X 1 Ai) 
Re (A) 9 (°) IP (eV) 



LDA 


1.03 


105.0 


6.2 


GW 






10.7 


Expt. [11, 12] 


1.01 


106.7 


10.1 






Excited state {A 1 A'J,) 






Re (A) 


o (°) 


T e (eV) 


CLDA 


1.08 


120 


5.05 


GW-BSE 


1.08 


120 


5.52 


CASSCF [13] 


1.06 


120 


5.49 


CEPA [13] 


1.06 


120 


5.63 


Expt. [12] 


1.08 


120 


5.7 



Table 2: Ground and excited state data for NH 3 : equilibrium bond length R e , H-N-H angle 6, ionization po- 
tential IP, and X x Ai — > A 1 A^ transition energy T e . The experimental T e contains a zero-point contribution 
of unknown but presumably small size. 
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